Availability and Quality Analysis Plot and Discussion -- Dataset for the SimAUD 2015 Conference
The goal of this analysis is to create plots to show the availability and quality of the data from the open dataset.
First, we'll look at the power meter data:
In [1]:
import pandas as pd
from matplotlib import ticker
In [2]:
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.backends.backend_pdf import PdfPages
import os
import matplotlib.dates as mdates
import datetime as dt
Load the raw datafile:
In [3]:
df = pd.read_csv("./SimAud_OpenDataSet.csv", parse_dates=True, index_col="Date Time")
Explore the dataset
In [5]:
df.info()
In [10]:
df.resample('W').plot(figsize=(20,300), subplots=True)
Out[10]:
Detect NaN's and "Flat lined" data
- Flat-lined data will be defined in this case as data in which the relative std dev is below a certain threshold
Let's look at a set of 10 streams as a example:
In [11]:
testcases = df.ix[:,1:10].resample("W")
In [12]:
testcases.plot(figsize=(20,10), subplots=True)
Out[12]:
Let's use the pandas tranform function to change the values for each day into the following ordinal values
- NaN -- 0
- 0 Values -- 1
- Outliers -- 2
- "Good Data" -- 3
First, filter out the outliers because they're dependent on the distributions
In [13]:
for datapoint in testcases.columns:
print testcases[datapoint][np.abs(testcases[datapoint]) > (testcases[datapoint].mean() + 3 * testcases[datapoint].std())]
In [14]:
outliers = pd.DataFrame()
for datapoint in testcases.columns:
# #The outliers first = 2
# testcases[datapoint][np.abs(testcases[datapoint]) > (testcases[datapoint].mean() + 3 * testcases[datapoint].std())]) = 2
# #Then non-outliers = 3
# print datapoint
tempset = pd.DataFrame(testcases[datapoint][np.abs(testcases[datapoint]) > (testcases[datapoint].mean() + 3 * testcases[datapoint].std())])
outliers = pd.merge(outliers, tempset, right_index=True, left_index=True, how='outer')
In [15]:
outliers[outliers.isnull() == False] = 2
Now, let's filter out the NaN, 0, and all other values
In [16]:
zeros = testcases[testcases == 0]
NaNs = testcases[testcases == np.nan]
In [17]:
zeros.info()
Transform the dataframe into the discrete variables
In [18]:
testcases[testcases>0].info()
In [19]:
testcases[testcases > 0] = 3
testcases[testcases <= 0] = 1
Update the dataframe with outliers extracted previously: http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.update.html
In [20]:
testcases.update(outliers)
testcases.fillna(value=0, inplace=True)
In [21]:
testcases.plot(figsize=(20,10), subplots=True)
Out[21]:
To get the discrete colors: http://gotoanswer.stanford.edu/matplotlib_discrete_colorbar-8699380/
Positioning colorbar: http://stackoverflow.com/questions/13310594/positioning-the-colorbar
Latexify: http://nipunbatra.github.io/2014/08/latexify/
Shift the plot a bit: http://stackoverflow.com/questions/6774086/why-is-my-xlabel-cut-off-in-my-matplotlib-plot
In [22]:
SPINE_COLOR = 'gray'
In [23]:
SPINE_COLOR = 'gray'
def latexify(fig_width=None, fig_height=None, columns=1):
"""Set up matplotlib's RC params for LaTeX plotting.
Call this before plotting a figure.
Parameters
----------
fig_width : float, optional, inches
fig_height : float, optional, inches
columns : {1, 2}
"""
# code adapted from http://www.scipy.org/Cookbook/Matplotlib/LaTeX_Examples
# Width and max height in inches for IEEE journals taken from
# computer.org/cms/Computer.org/Journal%20templates/transactions_art_guide.pdf
assert(columns in [1,2])
if fig_width is None:
fig_width = 3.39 if columns==1 else 6.9 # width in inches
if fig_height is None:
golden_mean = (sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_height = fig_width*golden_mean # height in inches
MAX_HEIGHT_INCHES = 8.0
if fig_height > MAX_HEIGHT_INCHES:
print("WARNING: fig_height too large:" + fig_height +
"so will reduce to" + MAX_HEIGHT_INCHES + "inches.")
fig_height = MAX_HEIGHT_INCHES
params = {'backend': 'ps',
'text.latex.preamble': ['\usepackage{gensymb}'],
'axes.labelsize': 10, # fontsize for x and y labels (was 10)
'axes.titlesize': 10,
'text.fontsize': 10, # was 10
'legend.fontsize': 8, # was 10
'xtick.labelsize': 8,
'ytick.labelsize': 10,
'text.usetex': True,
'figure.figsize': [fig_width,fig_height],
'font.family': 'serif'
}
matplotlib.rcParams.update(params)
def format_axes(ax):
for spine in ['top', 'right']:
ax.spines[spine].set_visible(False)
for spine in ['left', 'bottom']:
ax.spines[spine].set_color(SPINE_COLOR)
ax.spines[spine].set_linewidth(0.5)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
for axis in [ax.xaxis, ax.yaxis]:
axis.set_tick_params(direction='out', color=SPINE_COLOR)
return ax
In [24]:
latexify(fig_height=6)
numberofplots = 1
fig = plt.figure()
# Get the data
x = mdates.drange(testcases.index[0], testcases.index[-1] + datetime.timedelta(weeks=1), dt.timedelta(weeks=1))
y = np.linspace(1, len(testcases.columns), len(testcases.columns))
# Plot
ax = fig.add_subplot(numberofplots, 1, 1)
data = np.array(testcases.T)
cmap = cm.get_cmap('Reds', 4)
qmesh = ax.pcolormesh(x, y, data, cmap=cmap)
cbar = fig.colorbar(qmesh, ax=ax)
cbar.ax.tick_params(labelsize= 12)
tick_locator = ticker.MaxNLocator(nbins=3)
cbar.locator = tick_locator
cbar.update_ticks()
ax.axis('tight')
# try:
# plt.title(MasterPointList[(MasterPointList.Name == mappingname)].Name.iloc[0] + "/" +
# MasterPointList[(MasterPointList.Name == mappingname)].English_Description2.iloc[0] +
# "/"+ MasterPointList[(MasterPointList.Name == mappingname)].Description2.iloc[0], fontsize=26)
# except IndexError:
# continue
# Set up as dates
ax.xaxis_date()
#ax.yaxis.set_ticks()
fig.autofmt_xdate()
fig.subplots_adjust(hspace=.5)
format_axes(ax)
Out[24]:
Let's try it with all the points
In [25]:
df_daily = df.resample('W')
In [26]:
outliers = pd.DataFrame()
for datapoint in df_daily.columns:
tempser = pd.DataFrame(df_daily[datapoint][np.abs(df_daily[datapoint]) > (df_daily[datapoint].mean() + 3 * df_daily[datapoint].std())])
outliers = pd.merge(outliers, tempser, right_index=True, left_index=True, how='outer')
outliers[outliers.isnull() == False] = 2
df_daily[df_daily > 0] = 3
df_daily[df_daily <= 0] = 1
df_daily.update(outliers)
df_daily.fillna(value=0, inplace=True)
Output the results for replication without having to redo the above analysis
In [27]:
df_daily.to_csv('AvailabilityAnalysisResuldf_daily.csv')
Ok, now plot it
In [44]:
latexify(fig_height = 6)
numberofplots = 1
fig = plt.figure()
# Get the data
x = mdates.drange(df_daily.index[0], df_daily.index[-1] + datetime.timedelta(weeks=1), dt.timedelta(weeks=1))
y = np.linspace(1, len(df_daily.columns), len(df_daily.columns))
# Plot
ax = fig.add_subplot(numberofplots, 1, 1)
data = np.array(df_daily.T)
cmap = cm.get_cmap('Reds', 4)
qmesh = ax.pcolormesh(x, y, data, cmap=cmap)
cbaxes = fig.add_axes([0.2, 0.1, 0.7, 0.02])
cbar = fig.colorbar(qmesh, ax=ax, orientation='horizontal', cax = cbaxes)
cbar.ax.tick_params(length = 0)
tick_locator = ticker.MaxNLocator(nbins=3)
cbar.locator = tick_locator
cbar.update_ticks()
cbar.set_label('Quality Metric')
ax.axis('tight')
# try:
# plt.title(MasterPointList[(MasterPointList.Name == mappingname)].Name.iloc[0] + "/" +
# MasterPointList[(MasterPointList.Name == mappingname)].English_Description2.iloc[0] +
# "/"+ MasterPointList[(MasterPointList.Name == mappingname)].Description2.iloc[0], fontsize=26)
# except IndexError:
# continue
# Set up as dates
ax.xaxis_date()
#ax.yaxis.set_ticks()
fig.autofmt_xdate()
fig.subplots_adjust(hspace=.5)
ax.set_xlabel("Timeline")
ax.set_ylabel("Number of Points")
ax.set_title("kWh Meter Quality Screening")
plt.tight_layout()
plt.subplots_adjust(bottom=0.22)
format_axes(ax)
plt.savefig("availability_map.png", dpi=500)
Let's sort and then replot
In [29]:
df_daily_sorted = df_daily.T.ix[df_daily.T.sum(axis=1).order().index]
In [30]:
df_daily_sorted.head()
Out[30]:
In [31]:
latexify(fig_height = 6)
numberofplots = 1
fig = plt.figure()
# Get the data
x = mdates.drange(df_daily_sorted.T.index[0], df_daily_sorted.T.index[-1] + datetime.timedelta(weeks=1), dt.timedelta(weeks=1))
y = np.linspace(1, len(df_daily_sorted.T.columns), len(df_daily_sorted.T.columns))
# Plot
ax = fig.add_subplot(numberofplots, 1, 1)
data = np.array(df_daily_sorted)
cmap = cm.get_cmap('Reds', 4)
qmesh = ax.pcolormesh(x, y, data, cmap=cmap)
cbaxes = fig.add_axes([0.2, 0.1, 0.7, 0.02])
cbar = fig.colorbar(qmesh, ax=ax, orientation='horizontal', cax = cbaxes)
cbar.ax.tick_params(length = 0)
tick_locator = ticker.MaxNLocator(nbins=3)
cbar.locator = tick_locator
cbar.update_ticks()
cbar.set_label('Quality Metric')
ax.axis('tight')
# try:
# plt.title(MasterPointList[(MasterPointList.Name == mappingname)].Name.iloc[0] + "/" +
# MasterPointList[(MasterPointList.Name == mappingname)].English_Description2.iloc[0] +
# "/"+ MasterPointList[(MasterPointList.Name == mappingname)].Description2.iloc[0], fontsize=26)
# except IndexError:
# continue
# Set up as dates
ax.xaxis_date()
#ax.yaxis.set_ticks()
fig.autofmt_xdate()
fig.subplots_adjust(hspace=.5)
ax.set_xlabel("Timeline")
ax.set_ylabel("Number of Points")
ax.set_title("kWh Meter Quality Screening")
plt.tight_layout()
plt.subplots_adjust(bottom=0.22)
format_axes(ax)
plt.savefig("availability_map_sorted.png", dpi=500)
Cool, now let's zoom on a specific datapoint stream and show it individually with the colors
Let's find the best candidate by doing a subplot of all the points
In [32]:
df_daily.plot(figsize=(20,400), subplots=True)
Out[32]:
Pick a point - 1169
In [33]:
targetedpoint = df["1169"]
targetedpoint.resample('D').truncate(before='2011-12-01', after='2012-06-01').plot()
Out[33]:
In [34]:
targetedmetric = df_daily["1169"]
targetedmetric.truncate(before='2011-12-01', after='2012-06-01').plot()
Out[34]:
In [50]:
latexify(columns=1)
numberofplots = 1
fig = plt.figure()
# Plot
ax = fig.add_subplot(numberofplots, 1, 1)
y = targetedpoint.resample('D').truncate(before='2011-12-01', after='2012-06-01').values
x = np.arange(len(y))
t = targetedmetric.truncate(before='2011-12-01', after='2012-06-01').values
ax.plot(x, y, c='black', )
ymin, ymax = ax.get_ybound()
xmin, xmax = ax.get_xbound()
cmap = cm.get_cmap('Reds', 4)
im = ax.imshow(t.reshape(1, t.size), extent=[xmin,x.max(),ymin,ymax], alpha=.5, cmap=cmap, interpolation='none')
ax.set_xlabel("Timeline")
ax.set_ylabel("Point Number")
ax.set_title("kWh Meter Quality Screening")
plt.tight_layout()
plt.subplots_adjust(bottom=0.22)
format_axes(ax)
plt.savefig("availability_singleplot.png", dpi=1000)
plt.show()
different point
In [51]:
targetedpoint = df["806"]
targetedpoint.resample('D').truncate(before='2008-04-01', after='2009-02-01').plot()
Out[51]:
In [52]:
targetedmetric = df_daily["806"]
targetedmetric.truncate(before='2008-04-01', after='2009-02-01').plot()
Out[52]:
In [53]:
latexify()
numberofplots = 1
fig = plt.figure()
# Plot
ax = fig.add_subplot(numberofplots, 1, 1)
y = targetedpoint.resample('D').truncate(before='2008-04-01', after='2009-02-01').values
# x = mdates.drange(targetedpoint.index[0], targetedpoint.index[-1] + datetime.timedelta(days=1), dt.timedelta(days=1))
x = np.arange(len(y))
t = targetedmetric.truncate(before='2008-04-01', after='2009-02-01').values
ax.plot(x, y, c='black')
ymin, ymax = ax.get_ybound()
xmin, xmax = ax.get_xbound()
cmap = cm.get_cmap('Reds', 4)
im = ax.imshow(t.reshape(1, t.size), extent=[xmin,x.max(),ymin,ymax], alpha=.5, cmap=cmap, interpolation='nearest')
# ax.set_aspect(ax[0].get_aspect())
# ax.set_title(point)
# cax = fig.add_axes([0.2, 0.1, 0.7, 0.02])
# fig.colorbar(im, cax=cax, orientation='horizontal')
cbaxes = fig.add_axes([0.2, 0.2, 0.7, 0.02])
cbar = fig.colorbar(im, orientation='horizontal', cax = cbaxes)
cbar.ax.tick_params(length = 0)
tick_locator = ticker.MaxNLocator(nbins=3)
cbar.locator = tick_locator
cbar.update_ticks()
cbar.set_label('Quality Metric')
# # Set up as dates
# ax.xaxis_date()
# #ax.yaxis.set_ticks()
# fig.autofmt_xdate()
# fig.subplots_adjust(hspace=.5)
ax.set_xlabel("Days")
ax.set_ylabel("Energy [kWh]")
ax.set_title("Single Datastream Example")
plt.tight_layout()
plt.subplots_adjust(bottom=0.4)
format_axes(ax)
plt.savefig("availability_singleplot.pdf")
In []: